In [1]:
import sys
import os
from pathlib import Path

sys.path.append(Path.cwd().parent.joinpath("source").as_posix())

from utils import read_tiff, plot_sar, save_tiff, generate_graphs_from_pixel_coords
import plotly.offline as pyo
import plotly.graph_objects as go
import mvalab

# Set notebook mode to work in offline
pyo.init_notebook_mode()

import re
import numpy as np

from typing import Optional

import seaborn as sns
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams["figure.figsize"] = (20,20)

import random
random.seed(1)

from copy import deepcopy
from datetime import datetime

%load_ext autoreload
%autoreload 2
In [2]:
data_dir = "/mnt/hdd/SAR_images/processed/18MUU/cropped/descending"
image_paths = []
for file in sorted(os.listdir(data_dir)):
    if file.endswith(".tif"):
        filepath = os.path.join(data_dir, file)
        image_paths.append(filepath)

Temporal series¶

Average temporal backscatter profiles over the whole image after S1Tiling using $\sigma_0$¶

Assuming that in each pixel we have $\sigma_0$ i.e. the backscatter signal, we can plot its average value through time over the whole image. Note that this does not mean much, it would be better to compute this temporal series over a smaller patch of interested (e.g. an area where there was deforestation throughout the period vs an area where the forest was kept intact).

In [3]:
# Quick visualization to start
im = read_tiff(image_paths[0])
plot_sar(im)
In [4]:
im.shape
Out[4]:
(500, 600)
In [5]:
im.mean(), im.std()
Out[5]:
(0.27473092070629895, 0.17743026132446124)

Express $\gamma_0$ in dB¶

In [6]:
epsilon = 1e-7
In [7]:
temporal_means = generate_graphs_from_pixel_coords(
    data_dir = "/mnt/hdd/SAR_images/processed/18MUU/cropped/descending",
    title = 'Average Mean gamma_0 over all pixels per date',
    x_title = "Date",
    y_title = "gamma_0 (dB)",
    gamma0 = True,
    epsilon = epsilon,
    pixel_coords = None,
    offsets = None,
    return_results = True
)

Average temporal backscatter profiles over the whole image after RCR (10log10 being taken)¶

In [8]:
temporal_means_rcr = generate_graphs_from_pixel_coords(
    data_dir = "/mnt/hdd/SAR_images/processed/RCR/cropped/descending",
    title = 'Average RCR over all pixels per date',
    x_title = "Date",
    y_title = "RCR (dB)",
    gamma0 = False,
    epsilon = epsilon,
    pixel_coords = None,
    offsets = None,
    return_results = True
)

Mean $\gamma_0$ and RCR¶

In [9]:
pattern = r"20[0-9]{6}"
x = [re.findall(pattern, path)[0] for path in image_paths]
x = [datetime.strptime(date, "%Y%m%d") for date in x]

fig = go.Figure()
fig.add_trace(go.Scatter(x=x, y=temporal_means,
                    mode='lines+markers',
                    name='mean gamma_0'))
fig.add_trace(go.Scatter(x=x, y=temporal_means_rcr,
                    mode='lines+markers',
                    name='mean RCR'))
fig.update_layout(title='Average Mean RCR and Mean gamma_0 over all pixels per date',
                   xaxis_title='Date',
                   yaxis_title='RCR / gamma_0 (dB)')
fig.show()

Small visualization¶

In [10]:
index_shadows = read_tiff('/mnt/hdd/SAR_images/processed/shadows/cropped/descending/detected_shadows_image_indexes.tif')

colors = np.array([[0,0,0]] + sns.color_palette("rocket",n_colors=int(np.max(index_shadows) +1)))
color_index_shadows = colors[index_shadows.astype(int) + 1]
plt.imshow(color_index_shadows, interpolation='nearest')
plt.show()
plt.imsave( '/mnt/hdd/SAR_images/processed/shadows/cropped/descending/colored_index.png', color_index_shadows)

After outliers filtering with threshold = 16¶

In [11]:
index_shadows = read_tiff('/mnt/hdd/SAR_images/processed/shadows/cropped/descending/filtered_shadow_image_threshold_16.tif')

colors = np.array([[0,0,0]] + sns.color_palette("rocket",n_colors=int(np.max(index_shadows) +1)))
color_index_shadows = colors[index_shadows.astype(int) + 1]
plt.imshow(color_index_shadows, interpolation='nearest')
plt.show()
plt.imsave('/mnt/hdd/SAR_images/processed/shadows/cropped/descending/colored_filtered_shadow_image_threshold_16.png', color_index_shadows)

Focus on smaller patches where deforestation happened¶

In [12]:
patch_deforested1 = index_shadows[85:110, 320:340]
cropped_coords1 = [85, 110, 320, 340]
plt.imshow(patch_deforested1)
plt.show()
In [13]:
patch_deforested2 = index_shadows[75:160, 250:470]
cropped_coords2 = [75, 160, 250, 470]
plt.imshow(patch_deforested2)
plt.show()

Take a look at $\gamma_0$ and RCR on the first tiny crop¶

In [14]:
temporal_means_patch = generate_graphs_from_pixel_coords(
    data_dir = "/mnt/hdd/SAR_images/processed/18MUU/cropped/descending",
    title = 'Average Mean gamma_0 over all pixels per date, over a patch that has been deforested',
    x_title = "Date",
    y_title = "gamma_0 (dB)",
    gamma0 = True,
    epsilon = epsilon,
    pixel_coords = cropped_coords1,
    offsets = None,
    return_results = True
)
In [15]:
temporal_means_rcr_patch = generate_graphs_from_pixel_coords(
    data_dir = "/mnt/hdd/SAR_images/processed/RCR/cropped/descending",
    title = 'Average Mean RCR over all pixels per date, using gamma_0, over a patch that has been deforested',
    x_title = "Date",
    y_title = "RCR (dB)",
    gamma0 = False,
    epsilon = epsilon,
    pixel_coords = cropped_coords1,
    offsets = None,
    return_results = True
)

Look at the same values but for one pixel only¶

In [16]:
patch_deforested1[12:15,7:12]
Out[16]:
array([[ 1.,  1.,  4.,  5., -1.],
       [ 3.,  1.,  1.,  4.,  4.],
       [ 2.,  3.,  3.,  0.,  1.]])
In [17]:
offsets = [85, 320]
#pixel_coords = [10, 15, 8, 12]
pixel_coords = [12, 15, 7, 12]
In [18]:
temporal_means_pixel_def = generate_graphs_from_pixel_coords(
    data_dir = "/mnt/hdd/SAR_images/processed/18MUU/cropped/descending",
    title = 'Average Mean gamma_0 over all pixels per date, over a patch that has been deforested',
    x_title = "Date",
    y_title = "gamma_0 (dB)",
    gamma0 = True,
    epsilon = epsilon,
    pixel_coords = pixel_coords,
    offsets = offsets,
    return_results = True
)
In [19]:
temporal_means_rcr_pixel_def = generate_graphs_from_pixel_coords(
    data_dir = "/mnt/hdd/SAR_images/processed/RCR/cropped/descending",
    title = 'Average Mean RCR over all pixels per date, using gamma_0, over a patch that has been deforested',
    x_title = "Date",
    y_title = "RCR (dB)",
    gamma0 = False,
    epsilon = epsilon,
    pixel_coords = pixel_coords,
    offsets = offsets,
    return_results = True
)
  • temporal_means_rcr: contains the RCR for the given pixel where deforestation happened (note that the RCR can be seen as the derivative of the backscattered signal)

Find pixel where no deforestation happened¶

In [20]:
patch_deforested1[12:15,16:18]
Out[20]:
array([[-1., -1.],
       [-1., -1.],
       [-1., -1.]])
In [21]:
offsets = [140, 200]
pixel_coords = [12,15,16,18]
In [22]:
temporal_means_pixel_no_def = generate_graphs_from_pixel_coords(
    data_dir = "/mnt/hdd/SAR_images/processed/18MUU/cropped/descending",
    title = 'Average Mean gamma_0 over a patch where no deforestation happened',
    x_title = "Date",
    y_title = "gamma_0 (dB)",
    gamma0 = True,
    epsilon = epsilon,
    pixel_coords = pixel_coords,
    offsets = offsets,
    return_results = True
)
In [23]:
temporal_means_rcr_pixel_no_def = generate_graphs_from_pixel_coords(
    data_dir = "/mnt/hdd/SAR_images/processed/RCR/cropped/descending",
    title = 'Average Mean RCR using gamma_0, over a patch where no deforestation happened',
    x_title = "Date",
    y_title = "RCR (dB)",
    gamma0 = False,
    epsilon = epsilon,
    pixel_coords = pixel_coords,
    offsets = offsets,
    return_results = True
)

Take a look at a lodged area (not only the border as previously)¶

In [24]:
patch_deforested1[12:15,2:5]
Out[24]:
array([[-1., -1., -1.],
       [-1., -1., -1.],
       [-1., -1., -1.]])
In [25]:
offsets = [140, 200]
pixel_coords = [12, 15, 2, 5]
In [26]:
temporal_means_pixel_lodged = generate_graphs_from_pixel_coords(
    data_dir = "/mnt/hdd/SAR_images/processed/18MUU/cropped/descending",
    title = 'Average Mean gamma_0 over a lodged area',
    x_title = "Date",
    y_title = "gamma_0 (dB)",
    gamma0 = True,
    epsilon = epsilon,
    pixel_coords = pixel_coords,
    offsets = offsets,
    return_results = True
)
In [27]:
temporal_means_rcr_pixel_lodged = generate_graphs_from_pixel_coords(
    data_dir = "/mnt/hdd/SAR_images/processed/RCR/cropped/descending",
    title = 'Average Mean RCR using gamma_0 over a lodged area',
    x_title = "Date",
    y_title = "RCR (dB)",
    gamma0 = False,
    epsilon = epsilon,
    pixel_coords = pixel_coords,
    offsets = offsets,
    return_results = True
)

Visualize the three area in the same graph¶

  • shadow: temporal_means_pixel_def
  • lodged: temporal_means_pixel_lodged
  • forest: temporal_means_pixel_no_def
In [28]:
pattern = r"20[0-9]{6}"
x = [re.findall(pattern, path)[0] for path in image_paths]
x = [datetime.strptime(date, "%Y%m%d") for date in x]

fig = go.Figure()
fig.add_trace(go.Scatter(x=x, y=temporal_means_pixel_def,
                    mode='lines+markers',
                    name='shadow'))
fig.add_trace(go.Scatter(x=x, y=temporal_means_pixel_lodged,
                    mode='lines+markers',
                    name='lodged area'))
fig.add_trace(go.Scatter(x=x, y=temporal_means_pixel_no_def,
                    mode='lines+markers',
                    name='intact forest'))
fig.update_layout(title='Temporal backscatter proļ¬les of an area being logged, the corresponding shadow (on one of its edges when deforestation occurs) and an intact forest',
                   xaxis_title='Date',
                   yaxis_title='gamma_0 (dB)')
fig.show()